##################################################################

# Gautam Nair
# gautam.nair@yale.edu
# Misperceptions of Relative Affluence and Support for International Transfers
# Make Table 6: Heterogeneous Effects By Priors

##################################################################

# setwd("")

##################################################################
# loading packages
##################################################################

library(Hmisc)
library(foreign)
library(sandwich)
library(lmtest)
library(numDeriv)
library(stargazer)
library(ggplot2)
library(plyr)
library(gridExtra)
library(ri)
library(dplyr)
library(plyr)
library(scales)

##################################################################

rm(list=ls())

##################################################################
# Heterogeneous Treatment Effects by Priors
##################################################################

data.working <- readRDS("d_r_cleaned_analysis_dataset.Rda")

data.working$interaction <- data.working$group3*data.working$rank.median.over.incorrect

dep.var <- c("q13ab.3","q13ab.3.positive", "q8.toolittle", "q11.yes")

ind.group <- c("group3" , "rank.median.over.incorrect", "interaction")

reg.models <-vector("list", length(dep.var)) 
se.models <-vector("list", length(dep.var))
matrix <- vector("list", length(dep.var))

x=0 

for(i in 1:length(dep.var)){
	
	my.formula <-paste(dep.var[i],'~', paste(ind.group, collapse= ' + '))
	temp.model <- lm(my.formula, data=data.working)
	temp.se <- coeftest(temp.model, vcov=vcovHC(temp.model,type="HC2"))[,2]
	x = x+1
	se.models[[x]] <- temp.se
	reg.models[[x]] <- temp.model

}

varlabels <-       c("Information", "Biased Perception", "Information x Biased Perception" , "Constant" )
spectitle <- c("Heterogeneous Effects by Priors")
outputfile <- c("tf_t_06_priors_heterogeneity.txt")
depvarcaption <- c("Heterogeneous Effects by Priors")
cutoffs.star <- c(0.05, 0.01)
char.star = c("**", "***")
columnlabels <-c("Perc.", "+ve Donation (1/0)", "Increase Budget (1/0)", "Cut Protections (1/0)")

stargazer(reg.models[[1]], reg.models[[2]], reg.models[[3]], reg.models[[4]],
          se=list(se.models[[1]], se.models[[2]], se.models[[3]], se.models[[4]]),
         type="text",
         title= spectitle,
          out=outputfile,
          no.space=TRUE, 
          model.numbers= TRUE,
         notes.append=FALSE,
          align=FALSE,
         covariate.labels=varlabels,
          header= FALSE,
          font.size="small",
          model.names=TRUE,
          digits=2,
          omit.stat = c("rsq", "f","adj.rsq"),
          column.labels  = columnlabels,
          dep.var.caption  = "",
          dep.var.labels.include=FALSE,
         df=FALSE,
         omit.table.layout="n",
                  star.cutoffs = cutoffs.star,
     star.char = char.star
)   